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Abstract We introduce a method for computing prob¬ 
abilities for spontaneous activity and propagation fail¬ 
ure of the action potential in spatially extended, conduc¬ 
tance-based neuronal models subject to noise, based on 
statistical properties of the membrane potential. We 
compare different estimators with respect to the quality 
of detection, computational costs and robustness and 
propose the integral of the membrane potential along 
the axon as an appropriate estimator to detect both 
spontaneous activity and propagation failure. Perform¬ 
ing a model reduction we achieve a simplified analyti¬ 
cal expression based on the linearization at the resting 
potential (resp. the traveling action potential). This al¬ 
lows to approximate the probabilities for spontaneous 
activity and propagation failure in terms of (classical) 
hitting probabilities of one-dimensional linear stochas¬ 
tic differential equations. The quality of the approxi¬ 
mation with respect to the noise amplitude is discussed 
and illustrated with numerical results for the spatially 
extended Hodgkin-Huxley equations. Python simula¬ 
tion code is supplied on GitHub under the link https:// 
github.com/deristnochda/Hodgkin-Huxley-SPDE. 
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1 Introduction 


Noise is an inherent component of neural systems that 
accounts for various problems in information processing 
at all levels of the nervous system, see e. g. the review 


Faisal et al. (2008) for a detailed discussion. In partic¬ 


ular, channel noise has been identified as an important 
source of various types of variability in single neurons. 
Examples are the noise induced phenomena as observed 


Faisal & Laughlin (2007). The timing of action po¬ 


tentials can be highly sensitive with respect to fluctua¬ 
tions in the opening and closing of ion channels leading 


to jitter and stochastic interspike intervals (Horikawa 


19911. This effect becomes important in thin axons with 


diameter of less than 1/rm. Furthermore, there appear 
stochastic patterns in the grouping of action potentials, 
and action potentials can vanish due to noise interfer¬ 
ence or spontaneously emerge without apparent synap¬ 
tic input. 

When it comes to the mathematical modeling of the 
membrane potential in axons, in particular in thin ones, 
channel noise therefore has to be taken into account. 
For a discussion and comparison of the various types 
of adding noise to conductance-based neuronal models 
such as the classical Hodgkin-Huxley equations we refer 
to Goldwyn & Shea-Brown (2011). Concerning spatially 


extended models, in e. g. Tuckwell & dost (2010 2011); 


Tuckwell (20081 it has been shown that already simple 


additive noise, uncorrelated in space and time, accounts 
for a large range of variability in the action potential. 
That includes variability in the repetitive generation of 
action potentials, deletion of action potentials or prop¬ 
agation failures and spontaneously emerging action po¬ 
tentials or spontaneous activity. Because of this observa¬ 
tion, we restrict ourselves to such a simple model of the 
noise that as a byproduct reduces the computational 
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and analytical complexity. However, the proposed de¬ 
tection and estimation method can be applied to more 
complex models and e.g. even the full Markov chain 
dynamics of channel noise can be used. 

It is the purpose of this work to introduce a method 
to compute in a mathematical consistent way the prob¬ 
abilities of those last two events. This is done for gen¬ 
eral spatially extended neuronal models with additive 
noise, both numerically and theoretically, in terms of 
statistical quantities of the membrane potential. A suit¬ 
able statistical estimator for such kind of characteristics 
should have the following desired properties: It is au¬ 
tomatically evaluable to do Monte-Carlo simulations; 
it strictly separates the considered event from different 
ones; it is a low dimensional function of the observ¬ 
ables; it is relatively robust to stochastic perturbations 
and uncertainty in the observables. We compare differ¬ 
ent estimators with respect to the quality of detection, 
computational costs and robustness. In order to further 
reduce the computational costs and to obtain a simpler 
analytical description, we perform a consistent model 
reduction, with respect to these statistical quantities, 
to a one-dimensional linear stochastic differential equa¬ 
tion that allows to compute the desired characteristics 
without necessarily simulating the full system. 

The method is illustrated in a case study using the 
Hodgkin-Huxley equations (Hodgkin & Huxley 19521 
with two distinct parameter sets. With spatial diffusion, 
this is a system of partial differential equations that can 
serve as a model for the propagation of action potentials 
in the neuron’s axon. In particular, depending on the 
size of the stimulus there exist pulse-like solutions (ac¬ 
tion potentials) to these equations propagating along 
the spatial domain. Using these equations, we estimate 
the probabilities of spontaneous activity and propaga¬ 
tion failure. Although we only focus on these two ex¬ 
amples, the methods presented here can be used for a 
broader range of problems, in particular, similar model 
reductions can also be performed in order to compute 
time jitter and the variability in grouping patterns of 
action potentials. 

In our setting, we consider a simple spatial geometry 
of the axon that is a cylindrical shaped fiber. Thus the 
relevant spatial domain is an interval [0, L]. We pro¬ 
pose 4>{u) ■= Jq u(x) dx as an estimator for the de¬ 
tection of spontaneous activity and propagation fail¬ 
ures. Here, u is the space(-time)-dependent observable 
whose solution is pulse-formed. In the cases at hand, 
this will be the membrane potential, ^(u) is the area 
under the pulse considered as a graph with respect to 
the space variable that has the following properties: 
It is easy to extract automatically from the numeri¬ 
cal simulations; it signihcantly separates the number 


(Hodgkin & Huxley 1952 


sets. With spatial diffusion 


of observed pulses; it is a linear functional of only one 
observable; stochastic perturbations, in particular ad¬ 
ditive noise that is white in space (or of low corre¬ 
lation length) should cancel out through integration. 
The events of spontaneous activity and propagation fail¬ 
ure can both be defined as threshold crossings of the 
quantity ^(m) and therefore easily be estimated using 
a Monte-Carlo simulation. The results can be found in 
Section 3. In Section 4, we do a model reduction for 
this quantity, only assuming a reasonable local stabil¬ 
ity of the pulse and resting solutions. In particular, we 
deduce one-dimensional Ornstein-Uhlenbeck processes, 
that captures both probabilities in particular for small 
noise intensities remarkably well. 


2 Hodgkin-Huxley type equations 

In this article, we consider a spatially extended con¬ 
ductance based neuronal model with a simple one di¬ 
mensional domain (0, L) approximating the axon. This 
is most accurate in the case of a long axon, shaped as 
a cylinder with constant diameter. Our examples com¬ 
bine a Hodgkin-Huxley type model with diffusive spa¬ 
tial coupling to describe the evolution of the membrane 
potential u(t, x) in time and space by a system of partial 
differential equations involving the dimensionless potas¬ 
sium activation, sodium activation and sodium inacti¬ 
vation variables n(t, x), m{t, x) and h(t, x), respectively. 
This typically reads as 

CrndtU — ARi 

- - E^a) - 9 l{u- Al) 

dtn = aniu){l - n) - l3niu)n, (I) 

9(771 = amiu){l - m) - (3m{u)m, 
dth = ah{u){l -h) - I3h{u)h. 


Here, Cm is the membrane capacitance in /iF/cm^, d 
the axon diameter in cm, Ri the intracellular resistiv¬ 
ity in I7cm, (/k, 5Na, <?l the maximal potassium, sodium 
and leak conductance in mS/cm^. To further specify 
units, all times are in ms, voltages in mV and distances 
in cm. These standard parameters from the original 
work of Hodgkin & Huxley (1952) are used through¬ 
out: Ri = 34.5, Cm = 1, ffK = 36, gua = 120, gL = 0.3, 
Ek = —12, Axa = 115 and Al = 10. Note that the 
membrane potential is shifted by 65mV compared to 
the original values. In order to be in the regime of thin, 
unmyelinated axons, we choose a diameter oid = 0.5pm 
for all simulations and consider an axon length of L = 
1cm. 
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2.1 Two parameter sets for the (in)activation variables 

Equation § is missing the coefficients determining the 
evolution of the (in)activation variables. In the stan¬ 


dard model following Hodgkin & Huxley (1952) these 
are 


a„(u)= 


Oimiu) = 


10-u 


100 (e io“ — l) 
25-u 


10(e^'io'' 


- 1 )’ 


n(n) = —e- 


Mu) = 

Pm{u) =4e~^, 

Phiu) = 


30--U • 

e 10 - 1-1 


chain ion channel dynamics, we use this form of noise in 
our study. The reason is twofold: First, already such a 
kind of noise can qualitatively account for all of the phe¬ 


nomena observed in e.g. Faisal & Laughlin (20071 and 


second, it allows further analysis due to its simplicity. 
Mathematically speaking, current noise is realized as a 
two-parameter white noise 77 that is defined in terms 
of a cylindrical Wiener process W such that 77 = W. 
W = (W(t))t>o is a function valued process that can 
be formally represented by the infinite series 


^(t)(x) ■■= ^ enix)Mt), 


In the following, we refer to this model as (HH). A 
second model (HH) with a different behavior can be 
obtained by slight modification. Set 


dm(u) 


36 — M 

10 (e^ - 1 )’ 


Mu) 


1 


21.5-u 

e 10 


-f 1 


that amounts to a change in the sensitivity of the sodium 
(in)activation rates, and leave the rest unchanged. The 
result is a neuron much less sensitive to input, i. e. with 
a higher firing threshold. In the next section, models 
(HH) and (HH) are used to illustrate the phenomenon 
of spontaneous activity and propagation failure, respec¬ 
tively. 


2.2 A mathematical model 


Noisy perturbations of equation 0 can be realized as a 
stochastic partial differential equation (SPDE) on the 
Hilbert space (iJ, H-H) = L^(0, L) with inner product 
(•,•). The variables u{t),n{t),m{t) and h{t) are then 
function valued, thus we omit the x dependence in the 
notation. 

For the spatial diffusion, define the Laplace oper¬ 
ator Au '■= dfu supplemented with Neumann bound¬ 
ary conditions. We choose a sealed end a.i x = L, i. e. 
dxu{t, L) =0 for all t > 0 and model the input signal to 
the axon via an injected current in form of a rectangular 
pulse 


dxu{t,0) = - 


'KtP ' 


J(t) = 


J, t<T* 
0, t>T* 


( 2 ) 


Here, T* < 00 is the duration and J > 0 the amplitude 
of the signal. 

The question of how to add noise to equation 0 has 
been studied in the literature, see e. g. |Goldwyn &: Sheie] 
Brown (2011). Although it has been shown that cur¬ 


rent noise, i. e. uncorrelated additive noise in the volt¬ 
age variable, does not accurately approximate a Markov 


where (/3n(^))nGN is a family of iid real valued Brownian 
motions and 



is an orthonormal basis of H. For f,gGH one can 
calculate the covariance of this process as 


E[{W{t),f){W{s),g)] = {tAs){f,g), 


thus E[r]{t, x)ri{s, y)] = S(t — s)S(x — y), i. e. no correla¬ 
tion in either time nor space. Thus formally speaking, 
a cylindrical Wiener process is time-integrated space- 
time white noise. Equation ([^ then reads as 


du{t) = [XAu{t) + f{u{t),n{t), m(f), h(t))] dt 
+ adW{t), 

= an{u{t)) (1 - n{t)) - Pn{u{t))n{t), ( 3 ) 

= a^{u{t)) (1 - m(t)) - I3^{u{t))m{t), 

= ah{u{t)) (1 - hit)) - f3h{uit))hit). 


Together with suitable initial conditions, in our case 
the equilibrium values {u*, n*, m*, h*), being it* = 0 for 
(HH) and it* « —0.820 for (HH), as well as 


X = 


axiu*) 


Mu*) + Mu*)'' 


X = n,m, h, 


we refer to Sauer & Stannat (2014) for well-posedness 
of equation ([^. 


2.3 Linear stability of pulse and resting state 

If one injects an input above a certain threshold, the 
solution of equation 0 rapidly approaches a travel¬ 
ing pulse like solution. Denote X = (u,n,m,h)'^, then 
numerical simulations show that this traveling pulse is 
well-approximated by a solution of the form Xit,x) = 
Xix — ct) for a fixed reference profile X and pulse speed 
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c as long as the pulse did not reach the boundary. Let 
us call this solution X{t). 

Without any external input, the system Q remains 
in equilibrium if started there. Denote by X* this con¬ 
stant (in time and space) solution to the equations. 

The phenomena of interest in this work directly cor¬ 
respond to the stability properties of those two solu¬ 
tions X and X* . Although this has only been shown 
for general stochastic bistable equations, see e. g. [St an- 1 
), we assume a linear stability condition that 
possible to be extended to the higher dimen¬ 
sional Hodgkin-Huxley system. This linear stability as¬ 
sumption is then only used in Section 4 for a model re¬ 
duction. For convenience of notation, denote equation 
(j^ in the following abstract form 

dA:(t) = {AX(t) + F{X{t))^ dt + ad'W{t), (4) 

where A = {A, 0,... )^, W = {W, 0,... )^ and F is the 
appropriate nonlinear part of the drift. Also, denote by 
% = the state space of Q. Then, we assume 

the following geometrical condition of Lyapunov type 

([A-fVF(X*)]h,h)«<-/c*||h||^, (5) 

implying that the resting solution is locally exponen¬ 
tially attracting in "H, i. e. linearly stable. Moreover 

{[A+XF{X{r))]h,h)u < + ( 6 ) 

for all t S [TqjT], where dj^{t) = X{t). Here Tq is the 
time until X is in pulse form and T denotes the time, 
when the pulse has reached the boundary. The latter 
condition can be interpreted geometrically as follows: 
once it is formed, the traveling pulse solution is locally 
exponentially attracting in the subspace Ft '■= {h G 
% : {h,dj^(t))-H = 0} C H that is orthogonal to the 
direction of propagation. 


nat 


(2014 


should be 


2.4 Numerical method 


SPDE (j^ is a reaction diffusion equation coupled to 
a set of equations without spatial diffusion. Thus, the 
main issue from a numerical perspective is the simula¬ 
tion of equations of the form 


du{t) = [XAu{t) + f(t, dt + a dW (t) 

with Neumann boundary conditions as in §• The nu¬ 
merical method chosen for the integration of such a 
SPDE is a finite difference approximation in both space 
and time, see Sauer & Stannat (2015 2014) for details. 
For the space variable x we use an equidistant grid (xt) 
of size Z\a; = ^/n and replace the second derivative by 


its two-sided difference quotient. Boundary conditions 
are approximated up to second order, using the artifi¬ 
cial points X-i and The time variable t is dis¬ 

cretized to {tj ) using At = 1/m and a semi-implicit Eu¬ 
ler scheme. Approximating the variable u in the point 
{xi,tj) yields the following scheme. 


- AtFa 


"h ^ 

UN,j + l = UN,j + 3y|(2uAr_ij + i — 2uN,j + l) 

“h - 


- Atf{uN,j) + Ul/ 


for 1 < j < W — 1, where Jj is the discrete applied cur¬ 
rent and is a sequence of iid A/'(0,1)- 


distributed random variables. For details on conver¬ 


gence of this scheme and error rates we refer to Sauer 


& Stannat (20151. 


3 Reliability of signal transmission 

Let us first specify numerical parameters. We use N = 
500 gridpoints, i. e. Z\a; = 0.02, and At = 0.01 to sim¬ 
ulate the equations. Using the input of height J = 
0.001/rA and length T* = 0.5, in both models (HH) 
and (HH) a pulse is formed at the left boundary, trav¬ 
eling to the right, see Figure 



X 


Fig. 1 The time evolution of u using (HH) at ti = 10 (solid), 
t 2 = 20 (dashed) and = 30 (dotted) for the deterministic 
pulse (cr = 0) and one perturbed by noise (cr = 0.5). 


The problem at hand is how the presence of noise 
affects the generation and reliability of transmission of 
action potentials in the axon, similar to the studies by 


Faisal & Laughlin 

(2 

007 

1 for the Hodgkin-Huxley equa- 

tions and 

Tuckwell 

(2008 

) for the FitzHugh-Nagumo 
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equations. In particular, this section concerns two dis¬ 
tinct phenomena observed in these two studies. Faisal & 
Laughlin found that in the (HH) model action potential 
propagation is very secure, but in certain cases there 
spontaneously emerge action potentials somewhere along 
the axon due to the effect of noise {spontaneous activ¬ 
ity). This is illustrated in Figure]^ where an exemplary 
trajectory of such an event can be found. 



X 


Fig. 2 A realization of the event spontaneous activity is 
given by the solid, light gray trajectory. The three plots are 
the membrane potential u using (HH) at times ti = 13.5, 
t 2 = 14.5, ts = 16.5 from top to bottom. For comparison we 
include a trajectory, where there are only fluctuations around 
the resting potential (dashed, dark gray). For all of them, 
a = 0.372. 


On the other hand, Tuckwell observed that a pri¬ 
mary effect of noise on the action potential can be a 
breakdown of the pulse without any secondary phenom¬ 
ena such as spontaneous activity {propagation failure). 
An illustration is given in Figure comparing a fail¬ 
ure to a stable pulse. The equations Tuckwell used to 
model the neuron are, of course, different to the work 
by Faisal & Laughlin, however this discrepancy is not 
due to the choice of the neuron model but rather due to 
the choice of the particular parameter values describ¬ 
ing the model. These are directly linked to the stability 
of the traveling pulse and resting state. Inde^, slightly 
modifying sodium (in)activation in model (HH), we can 
observe occurences of propagation failure but no spon¬ 
taneous activity. In this work, (HH) is always used to 
study spontaneous activity and (HH) for the propaga¬ 
tion failures, since these are the prominent phenomena 
in the respective dynamical system. 

We aim to propose a simple statistical estimator 
that allows for detection of both spontaneous activity 



X 


Fig. 3 A realization of the event propagation failure is given 
by the solid, light gray trajectory. The three plots are the 
membrane potential u using (HH) at times = 14.5, t 2 = 16, 
ts = 18 from top to bottom. For comparison we include a 
trajectory, where no propagation failure occurs (dashed, dark 
gray). For all of them, cr = 0.504. 


and propagation failures. A first educated guess might 
suggest that checking for certain threshold crossings of 
the maximum height of the membrane potential, i. e. 
sup,jg(Q^^) u{x) > 0, is a good choice. Note that such a 
criterion has been used in Faisal & Laughlin (20071 to 
detect arrival times of action potentials. However, we 
suggest a different method using the following linear 
functional of the (shifted) membrane potential. 


(l>{u) := / u{x) — u* da;. 

Jo 

This describes the area below the pulse of the mem¬ 
brane potential shifted by the resting potential u*. Note 
that we can always change variables so that in the fol¬ 
lowing we assume w. 1. o. g. u* = 0. We choose the esti¬ 
mator <!> over any other pointwise criterion as e. g. the 
supremum for the following reasons. First, ^ is a linear 
functional of only one observable. Second, the action 
potential is not a point charge that propagates along 
the axon but it is rather spread out along some part 
of it that may reach up to a few cm in length. Thus, 
a global criterion as imposed by is more reasonable 
than a pointwise one. Moreover local fluctuations due to 
the noise should have a less pronounced effect. Third, <1> 
is not sensitive to fluctuations in the phase of the trav¬ 
eling pulse, which will be explained in the discussion 
section. 

Consider the deterministic solution (i. e. a = 0) u 
that is a traveling pulse and denote ^ := <P{u). As long 
as the pulse is formed, this quantity should stay more or 
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less constant. In the following, with abuse of notation 
we use <?(u) := <P{u)/^. Concerning the example paths 
in Figure 1^ and|^we can look at the corresponding time 
evolution of the area see Figure]^ 


4 
3 

> 8 . 2 
1 
0 

0 5 10 ti t2 ta 19 

3 
2 

> 8 . 1 

0 

-1 

0 5 10 ii t2 ta 20 

t 

Fig. 4 The evolution of the area (Top) For the same 
realizations as in Figure using (HH). Spontaneous activ¬ 
ity (light gray) and fluctuations around the resting potential 
(dark gray). (Bottom) For the same realizations as in Fig¬ 
ure [fusing (HH). Pulse with propagation failure (light gray), 
pulse without failure (dark gray). 





3.1 Spontaneous activity 

Since the estimator (p reliably discriminates between 
no, one or more pulses, it can be used to observe the 
probability of emerging secondary pulses. In this sce¬ 
nario, starting the model (HH) at the resting potential 
without any input signal through the Neumann bound¬ 
ary condition, we observe the solution for the time T 
the deterministic pulse u would need to reach the right 
boundary. For a given critical value 9 we dehne the 
event sup(g[o,T] ^ 9 as spontaneous activity 

for the noise amplitude a. Similar, the probability of 
spontaneous activity is 




sup <P(u'^{t)) > 9 
te[o.T] 


In this definition, the threshold 9 still has to be speci¬ 
fied. Experience with different parameter sets and other 
neuron models have shown that a suitable threshold de¬ 
pends heavily on these. Suitable is used here in the sense 
that the estimator indeed detects an emerging action 
potential when there is one. 

In the following we use T = 60 and M = 10 000 
realizations of (HH) to estimate s^. Figurej^shows that 
the curve u i—f So- shifts to the right as 9 is increased 
and stays unchanged for 9 > 0.52, which is in this case 
the suitable threshold to detect spontaneous activity. 



Fig. 5 Plot of So- vs. cr for different threshold values using 
the model (HH). 


3.2 Propagation failure 


Obviously, we can use <l> the other way rjamd to detect 
a propagation failure using the model (HH). Thus, we 
are in principle able to easily reproduce and general¬ 
ize the observations made in Tuckweli| (20081 in terms 
of variation of parameters, models and the number of 
Monte-Carlo realizations. Let Tg > 0 be a given, fixed 
initialization time until the pulse is formed. Also, recall 
that T denotes the time when the pulse has reached 
the boundary. Given a threshold 9 we define the event 
supjg[j.|j 7 .] ^(u‘^(t)) — $ > 9 as a propagation failure 
for the noise amplitude a. Similar, the probability of 
propagation failure is 


Pa ■ = 


sup <P{u'^{t)) - S > 9 

te[To,T] 


Remark 1 Numerically the stopping time T is imple¬ 
mented as follows. The axon is extended using a noise¬ 
less cable at the right boundary that allows to keep 
track of the pulse even if it already has left the original 
part of the axon. Applying the estimator <P on both the 
noisy and noiseless part makes it possible to determine 
whether and when a pulse has successfully reached the 
axon terminal. With this, we can compute a reference 
value to evaluate the quality of the estimator <P. 

With To = 10 Figure shows the probability of propa¬ 
gation failure Pa versus a for different threshold values 
compared to As 9 decreases, the curves converge 
to the reference curve. In particular, 9 = 0 seems like a 
suitable threshold in this scenario. 


4 Model reduction 

Obtaining an analytical expression for po- and Sa- is out 
of reach, considering these are the exit time probabili¬ 
ties of a nonlinear inhnite dimensional problem. How¬ 
ever, one can use the linear stability assumptions of 
both pulse and resting state made in Section 2.3 to ob¬ 
tain a simplified model. In this part we show that a 
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Fig. 6 Plot of Pct vs. a for different values of 6 using the 
model (HH). 


model reduction is indeed possible and propose a sim¬ 
ple, one-dimensional equation that mimics the behavior 
of the original problem and is able to capture the de¬ 
sired quantities, such as the probabilities of propagation 
failure and spontaneous activity. This has the following 
implications: First, the computational costs are reduced 
and second, we obtain a simplified analytical expression 
in terms of classical, known quantities. 

In view of assumption (§ our arguments for the use 
of ‘P can be strengthened by a simple observation. Let 
lu = (1,0,... )^ be the constant function equal to 1 in 
the M-component, then 

{lu,d^{t))'H = = 0 

for t G [To,T] since the integral is invariant to transla¬ 
tion of the pulse. Thus, 1„ € J-t for all t € [To,T]. 

The implications of this are the following. Consider 
the solution Z(t) to the linearization of @ neglect¬ 
ing all higher order terms. In particular, Z(t) is an 
Ornstein-Uhlenbeck process on "H. Writing T{t,s) = 
exp[J^* A + \7F{X{r)) dr] for the exponential of the lin¬ 
ear operator, the solution can be written using Duhamel’s 
principle as 

Z{t)=T{t,0)Z{0)+a [ T(t,s)dW(s). 

Jo 

Z{t) is a Gaussian process, uniquely characterized by 
its mean and variance 

E[z(t)] = r(t,o)z(o), 

Var[Z(t)] = j T{t, s)T{t, s)* ds, 

Jo 

where * denotes the adjoint operator. Now, recall P{u) = 
/g udx, hence P{u{t) — u{t)) = {u{t) — u{t),l)H ~ 
{Z{t), 1„)^. In particular, this is a linear functional of 
Z{t). Since Z is Gaussian, so is {Z(t), 1 u)'h with mean 
and variance 

E[{zit),iu)H] = (r(t,0)z(0),i,)«, 

Var[(Z(t),l„)„] =^2 f (r(t,s)T(t,s)*l„,l„)«ds. 

Jo 


Now, recall (§, i.e. the linear stability assumption for 
the pulse state. Note that G T^, i.e. orthogonal 
to the direction of pulse propagation, and therefore 
— 0. Hence, the linear operator T{t,s) 
satisfies the following inequality: 

||T(t,s)l,||«<e-«(‘-*)||l,||«. 

In particular it follows that 

E[{Z{t),K)n]<e--^\\zm\n\\lu\\n 
< VLe-^*\\zmn. 

Of course, this implies E[{Z{t), 1u)h] 0) which is one 

of the main advantages of choosing the estimator P. In 
contrast to this, the squared L^-norm \\u{t) — w(t)|||^ 
or also sup^jg^g j^) \u(t, x) — u{t, x)| might also serve as a 
measure of how close u is to the pulse solution. However, 
both will not converge to 0, since due to the noise u 
will never be adapted to the right phase of it. In our 
approach, we integrate the difference u — u with respect 
to a function orthogonal to the direction of propagation, 
hence our estimator does not perceive any phase shift 
and is locally exponentially stable around 0. Concerning 
the variance, we compute 

Var[(Z(t), 1„)«] = cr^ [ ||T(f,s)l„||^ds 

<a^ j\-^^^ds\\U\l<^. 

With the considerations above, the following Ansatz for 
a scalar valued stochastic differential equation for P is 
reasonable. 

dP(u(t) — u{t)) = —aP{u{t) — u{t)) dt -|- ad/3{t), 

where /3(t) := vT ^(IF(f),l) H defines a real-valued 
Brownian motion and a := y/La. Using linearity of P^ 
p := P[u(t)) and P{t) := P{u{t)) it follows that 

dP{t) = a(P — P{t)) dt + a dj3{t), P{Q) = P (7) 

is the approximating dynamics, a simple, one-dimensional 
Ornstein-Uhlenbeck process around the mean P. Also, 
Po- can be approximated by the exit time probability 


Pa ■ = 


sup P"’{t)-P>9 

te[To.Ti] 


Ti = E[T], that is a first passage time of the Ornstein- 
Uhlenbeck process. These are intensively studied in re¬ 
lation to stochastic LIF neurons, see Alili et al. (2005); 


Sacerdote & Giraudo (2013), and are in addition easily 


accessible numerically. 
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In this Ansatz, the whole complexity of the SPDE 
dynamics is reduced to the parameter a and the solu¬ 
tion to Q can be written down explicitly as 

<?(t) = ^ -h (<I>(0) - ^)e-“‘ + a f d/3(s). 

Jo 

Assuming the validity of this linear approximation, which 
will be true for small a, we can estimate a using mean 
and variance of ^(t). In particular, 


El(l>(t)]=S+(<P(0)-S)( 


Var[^>(t)] = E 


a / e 


- —2as 


— a(t—s) 


d/3(s) 


ds = 


La^ 


(l-e-2“‘). 


Hence, Yai[<P{t)] —>• i' 0 -^/ 2 a as t —>■ oo can be used to 
estimate a for large t, in our simulations t = 45, thus 
the difference to the limit is negligible. We apply the 
standard variance estimator 


VarM 


M M 

fe=i fc=i 


for cr = 0.024, the smallest a used in the simulations 
before. We arrive at 

La^ 


Q-m ■= 


20Var^ 


0.404, 


( 8 ) 


with again M = 10 000 realizations. 

Using the linearization around X* and the same 
Ansatz, we propose a similar Ornstein-Uhlenbeck pro¬ 
cess, whose hitting probabilities approximate s^. With 
<P(t) := (l>{u{t)) = {u{t), 1 )h and, of course, <P{u*) = 0 
this reads as 


d<?(t) = —/3^(t) dt + crd/?(t), <?(0) = 0. (9) 

Also, E[<?(t)] = 0 and Var[^(t)] = — and 

we estimate the rate fi using a = 0.012 via 


Pm 


La^ 


20Var' 


0.334 


( 10 ) 


M 


with M = 10 000 realizations. Figurej^shows the prob¬ 
abilities Pct and 


^<7 


sup > 9 

tG[To,Ti] 


as a function of cr for different thresholds 6 compared to 
the probabilities obtained using the SPDE. Note that 
the approximation becomes worse as 9 and a increase, 
which is expected since then the solution approaches 
the other equilibrium state and the linearization is not 
valid anymore. 




Fig. 7 Top plot: vs. a for two threshold values 6 in com¬ 

parison to the corresponding values for using the model 
(HH). Bottom plot: vs. cr for two threshold values 6 in 

comparison to the corresponding values for using (HH). 


5 Discussion 

In this article, we have introduced a method to compute 
probabilities for spontaneous activity and propagation 
failure in a consistent way with underlying spatially ex¬ 
tended, conductance-based neuronal models, based on 
certain statistical properties of the membrane poten¬ 
tial. Since the action potential in the neuron’s axon is 
not a point charge, but rather spread out in space, we 
advertise the use of a non-local criterion such as the one 
using (p. It may be interesting to find out how the axon’s 
length and diameter influence the quality of detection, 
since these are the relevant parameters concerning the 
width of an action potential. 

A further reduction in computational costs and a 
simplified analytical description can be achieved per¬ 
forming a model reduction with respect to the chosen 
estimator in a consistent way with the underlying 
spatially extended neuronal model. This is based on its 
linearization at the resting potential (resp. the traveling 
action potential) and allows to approximate the prob¬ 
abilities for spontaneous activity and propagation fail¬ 
ure in terms of (classical) hitting time probabilities of 
one-dimensional linear stochastic differential equations. 
Since the linearization is valid only locally, the approxi¬ 
mations Pa and Sa become worse for growing 9 and cr as 
shown in Figure For reasonable small 9 and a how¬ 
ever, the hitting probabilities of the one-dimensional 
stochastic differential equations are a solid approxima¬ 
tion to the full nonlinear, infinite dimensional SPDE. 
On the other hand. Fig. 7 also shows that the model 
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reduction can be used to find upper bounds for resp. 
Pct over a considerably larger range of a. 

In this study, we used the modified model (HH) to il¬ 
lustrate propagation failures. Although|Faisal fc Laugh-] 


lin (2007) found action potential propagation to be very 


secure with less than 1% failures, we have shown that 
little change in parameters produce a dynamical system 
with a totally different behavior. More precisely, rising 
slightly the sodium inactivation rate as in the modi¬ 
fied Hodgkin-Huxley system {HH) lowers excitability 
of the neuron and increases the probability of prop¬ 
agation failure. It may even become the predominant 
feature over spontaneous activation, similar to the case 
of the FitzHugh-Nagumo system, see Tuckwell (2008). 
It would be a interesting to see whether this computa¬ 
tional fact could be confirmed in experiments. 

As generalizations, we may incorporate more gen¬ 
eral noise, e. g. as suggested in |Goldwyn fc Shea-Brown] 
(2011) for the Hodgkin-Huxley model, and study how 
this affects the signal transmission. Note, that in the de¬ 
velopment of this study we have used e. g. conductance 
noise as presented in Linaro et al. (2011). This does not 
qualitatively change the behavior concerning pa- and Scr, 
but should be analyzed in comparison to the results in 
Faisal & Laughlin (2007| for the Hodgkin-Huxley equa¬ 


tions with ion channels modeled via Markov chains. Fu¬ 
ture work will also be concerned with the effect of noise 
on the generation of repetitive spiking, see Tuckwell & 


dost (2010), and the estimation of the speed of propa¬ 


gation. 


References 

Alili, L., Patie, P., & Pedersen, J. L. (2005). Rep¬ 
resentations of the first hitting time density of an 
Ornstein-Uhlenbeck process 1. Stochastic Models^ 
21 (4), 967-980. 

Faisal, A. A. & Laughlin, S. B. (2007). Stochastic sim¬ 
ulations on the reliability of action potential propa¬ 
gation in thin axons. PLoS Computational Biology, 
5(5), e79. 

Faisal, A. A., Selen, L. P., & Wolpert, D. M. (2008). 
Noise in the nervous system. Nature Reviews Neuro- 
seience, 5(4), 292-303. 

Goldwyn, J. H. & Shea-Brown, E. (2011). The what 
and where of adding channel noise to the Hodgkin- 
Huxley equations. PLoS Computational Biology, 
7(11), el002247. 

Hodgkin, A. L. & Huxley, A. F. (1952). A quantitative 
description of membrane current and its application 
to conduction and excitation in nerve. The Journal 
of Physiology, 777(4), 500-544. 


Horikawa, Y. (1991). Noise effects on spike propagation 
in the stochastic Hodgkin-Huxley models. Biological 
Cybernetics, 66(1), 19-25. 

Linaro, D. & Storace, M. & Giugliano, M. (2010). 
Accurate and fast simulation of channel noise in 
conductance-based model neurons by diffusion ap¬ 
proximation. PLoS Computational Biology, 7(3), 
el001102. 

Sacerdote, L. & Giraudo, M. T. (2013). Stochastic 
Integrate and Fire models: a review on mathemat¬ 
ical methods and their applications. In Stochastic 
Biomathematical Models (pp. 99-148). Springer. 

Sauer, M. & Stannat, W. (2014). Analysis and approx¬ 
imation of stochastic nerve axon equations. arXiv 
preprint arXiv:1402.4791, accepted for publication in 
Mathematics of Computation. 

Sauer, M. & Stannat, W. (2015). Lattice approxima¬ 
tion for stochastic reaction diffusion equations with 
one-sided lipschitz condition. Mathematics of Com¬ 
putation, 84{292), 743-766. 

Stannat, W. (2014). Stability of travelling waves 
in stochastic bistable reaction-diffusion equations. 
arXiv preprint arXiv:1404.3853. 

Tuckwell, H. G. (2008). Analytical and simulation 
results for the stochastic spatial Fitzhugh-Nagumo 
model neuron. Neural Computation, 20{12), 3003- 
3033. 

Tuckwell, H. G. & dost, J. (2010). Weak noise in neu¬ 
rons may powerfully inhibit the generation of repet¬ 
itive spiking but not its propagation. PLoS Compu¬ 
tational Biology, 6(5), el000794. 

Tuckwell, H. C. & dost, d. (2011). The effects of various 
spatial distributions of weak noise on rhythmic spik¬ 
ing. Journal of Computational Neuroscience, 30 {2), 
361-371. 




























